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1.  Introduction 

Consider  tbe  system  of  linear  equations 

Ax  ™  b  (1) 

where  the  coefficient  matrix  A  is  large  and  sparse  and  has  positive  definite  symmetric  part 
Af  wm  (A+At)f2.  In  this  paper,  we  compare  two  methods  for  solving  such  systems,  the 
generalised  conjugate  gradient  method  of  Concus  and  Golub  [2]  and  Widlund  [IOJ  and  Craig's 
method  (see  [0])  applied  to  a  symmetrically  preconditioned  auxiliary  system. 

Notation:  (y,z)  denotes  the  Euclidean  inner-product  y*t  and  ||-||  the  corresponding  noun.  If  Q 
is  a  symmetric  positive  definite  matrix,  then  (g,s)q  denotes  the  Q-inner  product  (Qy,z)  and  U*Hq 
the  corresponding  norm;  Q1^  denotes  any  square  root  of  Q;  and  denotes  (<j‘^*)  *.  Let 
A  —  M—N,  whence  —N  —  (A— A1) 1 2  is  the  skew-symmetric  part  of  A;  let  K  «■«  M~lN;  and  let 


""*■  "{(>+ 


m  —  0 
m  >  0 


2.  The  Generalised  Conjugate  Gradient  Method 

Concus  and  Golub  [2]  and  Widlund  (10]  proposed  tbe  Generalised  Conjugate  Gradient  (GCG) 
method  for  use  when  systems  of  the  form  Aft  —  4  are  “easy”  to  solve  (much  more  so  than  the 
original  system): 

Let  be  given  and  set  *(_i)  —  0. 

FOR  m  —  0  STEP  1  UNTIL  “CONVERGENCE”  DO 
rW  -  t  -  Ax{m) 

Solve  W”*  — 

.  -  ^»l) 


IcMislta  Par 

arts  mmi  ~ 
DTXC  TAB 
ifciaaaeoBcet 


7\ 


Distribution/ 


Availability  Codas 

fvall  and/or 
Special 


ft 


□  □ 


2 


-  *(*-»>  +  *m+i  |JM  +  *(-)  - 

The  cost  per  iteration  is  one  matrix  multiply  (by  .A),  one  solve  of  a  system  of  the  form  Mg  ■  4, 
and  2 n  multiplies. 

The  iterate  x^  can  be  characterised  as  the  unique  element  in  the  affine  Krylov  subspaee 

x{t)  +  Span{t/#),  Kv(*\  K*J*\ . K*_I«/,)}  »  *(#)+Sm 

satisfying  the  Galerkin  condition 


(*,  Ax(m)-h)  =  0  for  all  s  €  $m 
(see  [10]).  Moreover,  it  can  be  shown  that 
*("•>  =  *  +  xJKXx(t)-x)  , 

where  jrm(p)  w  even  (odd)  polynomial  of  degree  at  most  m  for  m  even  (odd)  and  *Jl) 
(see  [10]). 

The  iterate  can  also  be  characterised  as  the  best  approximation  from  a  certain 
m-dimensional  affine  subspace:1 


arg  min  Ijp-*#* 

»  €  *W+(/+Jf)Stt 

arg  min  Hr— 

,  e  ^+^+(i+K)su^ 


(see  [4]).  Equivalently, 


m  even  (=  2it) 
m  odd  (■»  2k+l) 


(2) 


where  Pm  is  the  space  of  all  real  polynomials  pjp)  of  degree  at  most  m  satisfying  pm(l)  —  1  and 
pm(-l)  -  (-1)-  (see  [4]).  Taking  p Jp)  -  TjiA~ln)/TJiA'1),  where  TJx)  is  the  mu 
Chebyshev  polynomial,  yields  the  error  bound* 


n**"0-*!*  < 


mm  +  i-m\~ 


where  R^A)  -  A'1  +  v/d_,+l  .  Taking  pjp)  -  #»*„_,(*)  yields  the  inequality 


1  However,  *<"*  m  nM  iht  boot  approximation  to  s  from  iho  ■natural*  affine  Krylov  Mbapnet  (sad  M) 

•  Tho  cm*  m  even  was  firwt  proved  by  flagman,  Luk,  sad  Young  M;  Widhrad  (!<|  givot  a  samawhM  washer 
hound. 


MM  **■'■**■ 
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A-i  <  A  for  all  m  >  1,  (4) 

which  shows  that  Me  even  and  odd  iterates  must  exhibit  the  tame  asymptotic  rate  of 
convergence  (cf.  [10]). 

3.  The  CSP  Method 

For  any  symmetric  positive  definite  matrix  Q,  the  system  (l)  is  equivalent  to  the 
symmetrically  preconditioned  system 

AS  s  [g-i/^g-i/*]  [g‘/*,]  -  [g-»/*|]  *  {  . 

If  we  apply  Craig’s  method  (see  [6])  to  this  auxiliary  system,  which  is  equivalent  to  applying  the 
conjugate  gradient  method  to  the  normal  equations 

AA'Q-b, 

then  the  resulting  method,  Craig’s  method  applied  to  the  Symmetrically  Preconditioned  auxiliary 
system  (CSP),  can  be  expressed  directly  in  terms  of  A,  z,  b,  and  Q  (see  [5]): 

LET  yW  (as  *<•>)  BE  GIVEN 

^ )  *  b  -  V#> 

Solve  Qr(,)  — 
pW-AVW 
Solve  Qpw~*pw 

FOR  k  —  0  STEP  1  UNTIL  “CONVERGENCE’’  DO 

tt  -  &  »a 

*  (p^^) 

^*+1) -V4)  +  q  *pw 
_  akApW 

Solve  Qt- w  «- 


Solve  Qp^l)  — 


„  (^‘)(  y^O) 

'*  (>s>V 
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The  cost  per  iteration  is  two  matrix  multiplies  (by  4  end  by  A1),  two  soiree  of  systems  of  the 
form  Mz  —  d,  and  5n  multiplies,  which  is  essentially  the  cost  of  two  GCG  iterations. 

The  iterate  can  be  characterised  as  the  unique  element  in  the  affine  Krylov  subspace 

*<•>  +  (g-‘4‘)  Span{F(,),  (Q-'AQ-'Ayi (Q-|A0“lA*f-|rw}  a  x(#)+Tm 

satisfying  the  orthogonality  condition 

(z,  ^~*)q  ■“  0  for  all  z  €  Tk. 

Thus 

If***  —  arg  min  ||y-:r|L  . 

»€*«+Tt 

and  the  standard  analysis  based  on  the  Chebysher  polynomials  yields  the  error  bound 

.  (S) 

/  +  p 

where  p  —  (ie+l)/(#e—  1)  and  k  is  the  condition  number  of  A. 

If  Q  —  M,  then 

*W+Tk  -  *<*>  +  (I+K)  Span {«/'*,  *V#>,  K**“V*)}  C  *«+(/+K)SM 

and 

“  arg  min  .  (®) 

»€»«+T* 

Moreover,  since*  it  ■■  t/l+d*,  the  error  bound  (5)  reduces  to  (3). 

But  an  even  stronger  relationship  exists  between  GCG  and  CSP.  If  m  is  even  (■  24),  then 
is  even  and  xm(l)  ■»  1  so  that 

»*(#»)  —  1  ~  (1+M)  Ptfc-|(#»)  (i—#*)  > 

where  n|t  is  an  even  polynomial  of  degree  at  most  24-2.  Thus,  by  equation  (2), 


*  tf  Um  somber  of  sqMtioM  is  im,  thta  y/i+A*  may  oafcr  bs  as  opper  bound  os  «. 


s 


*<**>  -  *(,)  -  (I+KbuJKXI-KX^-z) 

-xW  +  (I+K)ptk_t(Krf* 

e  *(,)+rt  ; 

i.e.,  the  beat  approximation  to  *  from  the  affine  svbspace  x^*4(/+liQ$i**  ^e#  'B  sm*^er 
affine  subspace  z(,4  7*.  But  since  is  the  best  approximation  to  *  from  *('474  (see  (6)),  it 
follows  that  x(tt)  —  i/A).  Hageman,  Luk,  and  Young  [8]  and  Elman  (5)  give  different  proofs  that 
the  two  methods  are  “virtually  equivalent.” 

The  cost  of  computing  ^  is  essentially  the  same  as  the  cost  of  computing  However, 
the  odd  iterates  generated  by  GCG  could  be  better  approximations  to  *  than  the  even  iterates 
(although  by  at  most  a  constant  factor  in  view  of  (4)).  Since,  in  addition,  GCG  requires 
somewhat  less  storage,  it  is  probably  the  better  method. 

4.  Two- Level  Methods 

But  what  if  systems  of  the  form  Ms  ■■  4  are  not  easy  to  solve?  Golub  and  Overton  [7)  have 
proposed  a  modification4  of  GCG  in  which  the  step 

SOLVE  Mi/",)  —  r*‘") 

is  replaced  by 

FIND  SOME  4m)  SATISFYING  |jWm)  -  » 

where  0  <  6  <  1  is  some  constant?  This  is  implemented  using  an  inner  iterative  method  to  find 
t/m)  on  the  mik  outer  iteration.  Basing  the  stopping  criterion  on  the  site  of  the  relative  residual 
has  the  effect  of  solving  Mt/"’  —  r1",)  to  increasing  absolute  accuracy  as  x(m)  converges  to  *. 

While  they  were  unable  to  analyte  this  two-level  scheme,  Golub  and  Overton  [7J  did  analyte 
a  similar  scheme  using  the  two-stage  Richardson  method  (also  a  three-term  recurrence)  as  the 
outer  iteration.  As  one  would  expect,  taking  i  closer  to  0  results  in  a  larger  number  of  inner 
iterations  per  outer  iteration  and  a  smaller  number  of  outer  iterations;  whereas  taking  f  closer  to 
1  results  in  a  smaller  number  of  inner  iterations  per  outer  iteration  but  a  larger  number  of  outer 
iterations.  The  same  behavior  for  the  two-level  GCG  method  can  be  seen  in  the  numerical 


*  Dsmbo,  Eieenstet,  sad  SteOuuc  [S]  stubs*  s  similar  modification  to  Newton's  method  ter  nonlinear  system*  of 
equations. 

1  The  csss  <  <■  0  corresponds  to  the  original  QCG  method. 


results  presented  in  Section  5. 

One  could  take  a  similar  approach  with  the  CSP  method.  Since  M  is  symmetric  and  positive 
definite,  a  logical  choice  for  the  inner  iteration  would  be  the  preconditioned  conjugate  gradient 
method  with  some  preconditioning  matrix  Q  (see  (1]).  Bnt  why  use  a  two-level  iteration  at  all 
when  one  can  simply  take  Q  «  Q  (instead  of  Q  Afft  The  numerical  results  presented  in 
Section  5  suggest  that  this  approach  is  superior. 

6.  Numerical  Results 

In  this  section,  we  reproduce  the  numerical  experiments  repented  by  Golub  and  Overton  [7] 
for  the  two-level  GCG  method  and  present  the  corresponding  results  for  the  CSP  method. 

Consider  the  elliptic  partial  differential  equation 
-An  +  (oo)^  +  auM  +  (tu)9  +  tuf  +  eu  —  / 
subject  to  Diriehlet  boundary  conditions  on  the  unit  square  fD,l]X[0,l],  where 
«(*,„)  -  5e*‘^,  *(*,»)  -  lOe*^ 

and  J[z,y)  and  the  boundary  conditions  are  chosen  to  make  the  solution  «(*,»)  ■■ 

The  five-point  centered  finite-difference  discretisation  on  a  rectilinear  grid  with  n  interior 
mesh  points  in  each  direction  leads  to  a  system  of  n*  linear  equations  in  which  Mu  corresponds  to 
— Au  +  eu.  Thus  we  use  the  fast  Poisson  solver  HWSCRT  from  FISHPACK  (9)  as  a 
preconditioning  for  CSP  and  for  an  inner  preconditioned  conjugate  gradient  iteration  in  the  two- 
level  GCG  scheme.  In  each  case,  the  stopping  criterion  was  <  1(T*. 

The  numbers  of  (outer)  iterations  and  Poisson  solves  are  given  in  Table  I  and  the  number  of 
Poisson  solves  is  plotted  against  6  in  Figure  1.  Clearly  CSP  is  a  better  method  than  GCG  for 
this  problem,  even  with  the  optimal  choice  of  6. 

I.  Conclusions 

/ 

As  we  have  seen,  if  systems  of  the  form  M*  "fare  “easy*  to  solve,  then  GCG  is  better  than 
CSP.  If  not,  then  CSP  is  superior.  Of  course,  it  is  not  clear  that  either  method  is  the  best 
possible  for  this  clam  of  problems. 
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Table  It  Number  of  (Outer)  Iterations  and  Poisson  Solves 

u  *  15  ft  s  81 
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a 

*P 

a 

.0 

33 

36 

.00100 

36 

266 

39 

279 

.00158 

36 

250 

39 

270 

.00251 

36 

247 

39 

261 

.00398 

37 

238 

41 

255 

.00631 

38 

226 

41 

244 

.0100 

38 

218 

42 

235 

.0158 

39 

209 

42 

216 

.0251 

40 

198 

43 

209 

.0398 

41 

193 

45 

196 

.0631 

44 

175 

47 

181 

.100 

52 

182 

71 

237 

.158 

135 

398 

273 

761 

89 

78 

42 

84 

a  =  nuebsr  of  (outsr)  iterations 
*.  =  nuebsr  of  Poisson  solves 
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